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1.  INTRODUCTION  AND  SUMMARY 


1.1  BACKGROUND 

DNA  (formerly  DASA)  became  involved  with  underground  nuclear  testing  in 
the  late  1960's  when  it  became  apparent  that  such  tests  could  be  useful  in  studying 
weapons  effects.  As  the  program  grew,  it  became  clear  that  planning  personnel  in 
the  field  needed  criteria  on  which  to  select  sites  for  future  tunnel  development. 
Earlier  site  separation  criteria  which  were  established^  were  chosen  in  a conservative 
way  because  the  data  available  at  the  time  were  limited.  However,  the  referenced 
criteria  were  adopted  as  the  standard  for  lack  of  any  rational  basis  to  change  them. 

As  a result  of  the  containment  program  supported  by  DNA  to  help  evaluate 
the  safety  of  their  underground  nuclear  tests,  a large  base  of  data  in  both  the 
material  properties  area,  as  well  as  in  the  area  of  reentry  observations,  has  been 
established.  In  addition  to  the  data  bases,  analytical  and  numerical  tools  for 
considering  various  aspects  of  containment  phenomenology  have  evolved. 

The  importance  of  siting  horizontal  line-of-site  nuclear  tests  in  materials 
which  are  at  once  low  in  shock  attenuation,  and  moderate  in  strength,  has  led  to 
extensive  material  properties  determinations  for  site  selection  purposes.  Some  areas 
of  Rainier  and  Aqueduct  mesas  are  clearly  better  than  others  according  to  these 
criteria,  and  in  order  to  maximize  the  use  of  these  areas  (which  are  considered  best 
for  containment)  a conflict  regarding  the  proximity  of  a proposed  test  to  an  old 
chimney  arises.  The  computations  reported  here  represent  a first  cut  at  understanding 
the  interactions  between  an  underground  test  and  existing  nuclear  explosion-produced 
chimney. 

1.2  STUDIES  CONDUCTED 

The  studies  which  have  been  conducted  are  of  two  general  types:  detonations 
adjacent  to  nearby  chimneys  and  detonations  at  the  working  point  of  a prior  explosion 
which  has  produced  a chimney.  The  nearby  chimney  studies  have  modeled  the 
chimneys  as  disk-shaped  regions  of  modified  materials  with  centers  located  at 
distances  of  122  meters  (400  feet)  and  152  meters  (500  feet)  from  the  working  point 
of  the  new  detonation.  For  each  location,  two  companion  computations  at  a nominal 
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10  kT  yield  were  performed  in  order  to  bound  the  effects  of  material  strength.  In 
all  calculations  the  volume  from  the  chimnied  cavity  was  spread  over  the  chimney 
volume.  The  strength  of  the  chimney  material  was  considered  to  be  the  same  as 
the  undisturbed  material  in  one  case  and  zero  (as  an  extreme  bound)  in  the  other. 
Preliminary  discussions  of  these  studies  have  already  been  published  as  References 
(2)  and  (3). 

The  compuations  performed  for  the  cases  where  detonation  is  to  take  place 
at  the  working  point  of  a prior  event  consider  the  effects  of  different  material 
properties  and  chimney  shape.  These  latter  studies  are  presented  for  the  first  time 
in  this  report.  A cylindrical  and  a conoidal  chimney,  each  presumed  to  have  been 
caused  by  10  kT  explosions  and  having  different  material  properties,  provide  the 
media  for  5 kT  explosions  placed  at  the  original  detonation  point.  The  properties 
of  the  cylindrical  chimney  material  include  a void  volume  from  the  original  burst 
cavity  distributed  through  the  chimney.  The  material  in  the  conoidal  chimney  was 
assumed  to  have  the  same  properties  as  the  undisturbed  tuffs  which  would  collapse 
into  chimneys  in  Rainier  Mesa.  This  material  has  a much  lower  air  void  content 
and  exhibits  higher  strength  than  that  in  the  cylindrical  chimney. 

1.3  SUMMARY  OF  RESULTS 

1.3.1  Pancake  Simulation 

The  disk-shaped  chimney,  "pancake",  centered  at  a range  of  152  meters  was 
found  to  have  very  little  influence  on  the  stress  fields  around  the  cavity  at  anytime, 
however,  a modest  decrease  in  axial  stress  compared  to  that  of  the  free-field  did 
occur  at  about  120  meters  near  the  immediate  vicinity  of  the  pancake  edge.  For 
the  pancake  centered  at  122  meters,  a 50  percent  decrease  in  axial  stress  was 
observed  immediately  in  front  of  the  edge  of  the  pancake  at  about  60  ms,  but  the 
cavity  region  experienced  only  about  a 15  percent  decrease  in  axial  stress  at  a time 
of  80  ms.  Neither  chimney  location  produced  more  than  a very  slight  change  in 
the  residual  stress  fields  around  the  cavity.  As  one  would  expect,  the  stress  field 
inside  both  chimneys  was  significantly  altered  from  that  of  the  free-field  due  to  the 
material  properties  differences. 

The  "fluid  filled"  chimneys  produced  a much  larger  effect  on  the  dynamic  and 


residual  stress  fields  than  did  the  chimneys  which  were  given  the  same  shear  strength 
as  the  surrounding  tuff.  Although  the  compressive  hoop  stress  was  only  slightly 
reduced,  the  thickness  of  the  residual  stress  region  around  the  cavity  was  modified 
by  the  presence  of  the  pancakes.  In  the  case  in  the  chimney  centered  at  122  meters, 
the  residual  hoop  stress  decreased  rapidly  at  about  65  meters  to  a "fractured"  region 
whereas  for  the  pancake  centered  at  152  meters,  the  rapid  drop  to  a fractured 
region  occurs  at  about  85  meters.  The  peak  residual  axial  stress  was  reduced 
approximately  15  percent  from  the  free-field  value  in  the  direction  of  the  far  chimney 
and  30  percent  towards  the  near  chimney. 

The  series  of  four  calculations  discussed  here  demonstrate  the  sensitivity  of 
the  dynamic  and  residual  stress  fields  to  the  material  properties  of  a nearby  chimney. 
The  pancake  geometry  is  inherently  conservative  regarding  the  propagation  of  reflec- 
ted waves  from  a nearby  chimney  back  towards  the  detonation  point.  Data  which 
has  become  available  subsequent  to  the  performance  of  these  calculations  on  the 
material  properties  of  fractured  and  reconstituted  tuff  samples  indicates  a strength 
reduction  of  approximately  25  percent  at  confining  pressures  above  1 kbar  and  30 
percent  at  0.5  kbar.  The  unconfined  strength  for  the  disturbed  material  was  measured 
to  be  zero,  but  the  strength  rose  rapidly  with  small  increases  in  confining  pressure. 
On  the  basis  of  these  test  results  and  the  recognition  that  at  burst  elevation  the 
chimney  tuff  is  recompacted  to  near  its  original  properties,  it  is  felt  that  the  "strong" 
chimney  material  gives  a better  estimate  of  the  chimney  response  than  the  zero 
strength  material. 

1.3.2  In-Chimney  Detonation 

The  most  striking  result  from  all  of  the  calculations  reported  here  was  the 
late  time  velocity  retained  by  the  chimney  material  for  the  cylindrical  chimney  case 
which  was  still  moving  at  350  ms.  For  the  5 kT  yield  used  for  the  in-chimney 
calculations  one  expects  the  free-field  stresses  to  settle  down  to  their  final  configu- 
ration at  about  200  ms.  While  the  hemispherical  region  below  the  cavity  was 
relatively  unperturbed  by  the  presence  of  the  chimney,  the  material  above  the  cavity 
had  hoop  stresses  below  that  of  the  magnitude  of  the  cavity  pressure  at  350  ms. 
Also,  a tensile  region  was  found  to  exist  at  that  time  which  formed  a valley  around 
the  region  of  high  residual  stress  and  intersected  the  tun  ?1  horizon. 


In  the  computation  of  the  conoidal  chimney  which  utilized  more  realistic 
material  properties  at  least  for  the  lower  region  of  the  chimney,  the  results  showed 
that  the  presence  of  the  chimney  created  small  perturbations  in  the  propagation  of 
the  shock  wave,  but  that  the  residual  stress  fields  at  a time  of  250  ms  were  only 
slightly  perturbed  in  the  direction  above  the  cavity.  While  a small  tensile  region 
at  the  edge  of  the  chimney  existed,  it  was  isolated  from  the  cavity  by  high  residual 
hoop  stresses. 

1.3.3  Conclusions 

Taken  together,  the  computations  for  both  the  pancake  and  the  in-chimney 
detonation  studies  show  the  importance  of  the  knowledge  of  material  properties  in 
chimney  regions.  For  the  pancake  simulations,  a chimney  centered  at  152  meters 
does  not  appear  to  create  significant  perturbations  to  either  the  dynamic  or  late-time 
stress  fields  from  an  adjacent  explosion  regardless  of  its  material  response.  Con- 
sidering the  conservative  assumptions  regarding  the  geometry,  it  is  expected  that  a 
chimney  centered  122  meters  from  an  explosion  of  the  magnitude  studied  here  would 
not  produce  significant  effects  unless  the  chimney  material  were  of  very  low  strength. 
In  order  to  study  the  influence  of  a nearby  chimney  of  very  low  strength,  the  simple 
pancake  geometry  would  have  to  be  abandoned  in  order  to  obtain  realistic  results. 

The  sensitivity  of  the  results  of  the  in-chimney  detonation  to  the  choice  of 
material  properties  suggests  that  an  extensive  data  acquisition  program  would  be 
required  to  justify  the  utilization  of  an  old  chimney  for  execution  of  another  event 
having  a yield  of  factor  of  2 lower  than  the  event  which  produced  the  chimney. 


2.  TECHNICAL  DISCUSSION 


2.1  OBJECTIVES 

The  objectives  of  the  work  reported  here  were  to  determine  whether  any 
dramatic  effects  which  could  influence  containment  would  occur  due  to  ground  shock 
interaction  with  prior  event-produced  chimneys  at  different  locations.  The  information 
available  regarding  the  strength  of  chimney  material  was  nonexistent  at  the  time 
the  studies  commenced.  However,  recognition  of  the  importance  of  strength  effects 
led  to  the  zero  strength  bounding  calculations  which  are  reported  here.  What  was 
anticipated  was  the  possible  occurrence  of  early-time  dynamic  effects  which  could 
cause  large  displacements  in  the  region  between  the  detonation  point  and  the  old 
chimney;  also  anticipated  were  the  possible  reduction  of  the  residual  stress  fields 
which  ordinarily  accompany  the  detonation  of  a nuclear  device  underground.  These 
stress  fields  are  considered  to  aid  in  containing  the  radioactive  gases  produced  by 
the  detonation. 

2.2  MODELING  THE  REAL  WORLD 

The  configurations  which  best  represent  situations  close  to  reality  are  generally 
three-dimensional  in  nature.  This  is  especially  true  when  one  considers  the  interaction 
of  an  underground  explosion  with  a chimney  produced  by  an  explosion  at  the  same 
elevation.  Furthermore,  geologic  layering  and  the  effects  of  gravity  are  additional 
aspects  of  the  phenomenology  which  if  included  would  dictate  that  three-dimensional 
models  would  be  required  for  analysis.  The  large  cost  of  exercising  three-dimensional 
computer  programs  is  well-known  and  the  hardware  available  generally  puts  severe 
limitations  on  the  resolution  with  which  one  can  compute  a given  region  of  space. 
For  these  reasons,  and  because  the  present  work  was  only  of  an  exploratory  nature, 
the  utilization  of  two-dimensional  axisymmetric  techniques  was  considered  appro- 
priate. The  STAR  code  which  was  used  for  all  the  computations  is  discussed  below. 

2.3  THE  STAR  CODE 

STAR  solves  the  finite  difference  equations  of  motion  of  compressible,  inviscid, 
continuous  media  in  a two-dimensional  Lagrangiari  framework.  This  approach  is 
particularly  advantageous  when  the  constitutive  relation  is  dependent  on  the  defor- 
mation history. 


Special  features  of  STAR  include  generalized  boundary  conditions,  much  faster 
energy  and  stress  iterations  than  in  early  2D  Lagrangian  codes, ^ an  improved 
artificial  viscosity  treatment,  an  improved  capability  to  detect  regions  requiring 
rezone,  an  automatic  rezoning  routine,  and  subcycling  of  grid  subregions.  In  addition 
to  the  person-time  saved  by  not  requiring  intervention  and  manual  operations  during 
runs,  automatic  rezoning  and  subcycling  techniques  can  save  factors  of  between  10 
and  100  in  computing  time  for  many  classes  of  problems. 

2.4  MATERIAL  RESPONSE  MODEL 

2.4.1  Equation  of  State  Library 

All  aspects  of  the  material  response  model  are  computed  in  a separate  library 
of  equation  of  state  routines  which  are  called  by  either  STAR  or  SIMONE.  This 
procedure  insures  that  identical  descriptions  of  the  material  behavior  are  employed 
in  comparable  one  and  two-dimensional  calculations.  This  equation  of  state  library 
permits  very  general  formulations  of  the  constitutive  response  of  real  materials 
including  the  response  of  porous,  visco-elastic  and  visco-plastic  materials.  Multi- 
material and  multi-equation  of  state  problems  are  readily  treated.  Energy  and 
pressure  dependent  yield  surfaces,  as  well  as  variable  elastic  moduli,  are  generally 
employed  for  most  materials  of  interest.  Although  not  used  here,  associative  flow 
rule  and  strain  dependent  yield  surface  formulations  are  available. 

2.4.2  Equation  of  State 

The  current  equation  of  state  library  contains  six  different  formulations  which 

specify  pressure  as  a function  of  density  and  internal  energy:  gamma-law  gas, 

Mie-Gruneisen,  polytropic,  Tillotson,  JWL  explosive,  and  a very  detailed  tabular 

(5) 

equation  of  state  for  water.  In  addition  to  these  response  formulations,  multi- 
material, tabular  equations  of  state  can  be  generated  by  combining  up  to  four 
materials  in  any  desired  ratio.  Any  of  the  six  basic  equations  of  state  can  be  used 
to  represent  the  response  of  each  component  of  the  mixture.  The  tabular  "isponse 
of  the  mixture  as  a whole  is  then  generated  automatically  by  iteration  to  achieve 
pressure  equilibrium  between  the  various  components.  Common  table  dimensions  are 
128  by  32  entries  in  u = (p/pQ-l)and  log(e)  space,  respectively,  spanning  a pressure 
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range  of  0.1  to  1 x 10°  MPa.  There  is  no  restriction  in  the  code  to  limit  the 

number  of  materials  using  the  multi-component  mixture  formulation  (other  than 

increased  storage  requirements)  and  calculations  with  multi-tabular  materials  have 
(6) 

been  made. 

In  both  the  chimney  and  pancake  studies  cavity  gas  was  modeled  as  a uniform, 

gamma-law  gas  with  y = 1.333.  For  the  regions  and  times  of  interest  in  these 

studies  it  is  well  known  that  detailed  calculations  of  the  device  interaction  with  the 

zero  room  and  surrounding  medium  are  not  required  to  satisfactorily  reproduce  the 
(7) 

ground  shock.  Additionally,  parametric  studies  which  model  cavity  gas  with  a 

chemical  equilibrium  equation  of  state  for  the  rock  vapor  suggested  that  the  gamma 

law  source  model  is  equivalent  to  a modest  increase  in  device  yield  in  comparison 

to  the  chemical  equilibrium  source,  and  that  effects  due  to  differences  in  the  timing 

of  the  source  pressure  decrease  are  of  second  order  by  comparison  to  the  effective 
(6) 

yield  reduction. 

2.4.3  Crush  Response 

Other  constitutive  models  in  addition  to  the  equations  of  state  discussed  in 
Section  2.4.2  are  required  to  treat  the  general  constitutive  response  of  porous 
materials.  For  such  materials  the  equation  of  state  is  assumed  to  described  only 
the  response  of  material  in  a fully  compacted  or  zero  void  condition.  The  response 
of  partially  crushed  material  is  then  referenced  to  the  fully  crushed  state  through 
the  distention  ratio, a,  defined  as 


ot  = v/v?  >1.0 


where  v and  vg  are  the  specific  volumes  of  the  porous  and  crushed  material, 
respectively.  The  pressure  of  the  porous  material  is  then  given  by 


p = - p (v  ,e) 
y a ys  ' s 


v-v 


(8) 

where  this  p-a  formulation'  ' is  completed  by  specifying  the  variation  of  <*  with  pg. 

The  functional  dependence  of  the  distention  ratio,  ct  , on  pressure  during 
loadir^j  can  be  derived  from  measured  p-v  response  data  provided  that  the  response 
of  a fully  compacted  sample  is  also  determined  over  the  same  pressure  range.  The 
resulting  crush  curve  has  frequently  been  approximated  in  earth  motion  codes  using 
a convenient  function  with  several  adjustable  coefficients  which  are  varied  to  give 
a "best  fit".  In  our  studies  the  crush  curve  is  represented  in  tabular  form  with  a 
a series  of  polynomial  fits  to  accurately  define  the  curve  between  tabulated  values. 
Thus  any  reasonable  crushup  behavior  can  be  modeled. 

In  principal,  the  unloading  response  can  also  be  measured  in  the  laboratory  . 
for  a family  of  unloading  paths  from  a series  of  points  on  the  crush  curve.  In 
practice,  this  information  is  seldom  available.  As  a consequence,  simple  models  for 
the  unloading  response  of  a have  generally  been  employed.  In  our  studies  release 
paths  from  pressures  above  the  elastic  limit,  Pg,  and  below  the  crush  pressure,  Pc, 
were  assumed  to  follow  paths  of  constant  a.  Load  paths  below  Pg  were  reversible, 
and  material  loaded  above  Pc  was  constrained  to  follow  the  solid  material  response 
with  a equal  to  the  fully  crushed  value  of  1.0. 

2.4.4  Elastic  Response 

At  low  stress  levels  competent  rock  media  behave  nearly  elastically.  In  this 
stress  regime  compressive  strain  paths  are  essentially  reversible,  so  that  neither 
dilatation,  nor  permanent  compaction  are  experienced  by  the  material,  and  its  response 
is  fully  described  by  the  various  elastic  moduli.  Thus  inclusion  of  this  low  stress 
behavior  in  the  numerical  model  is  relatively  straightforward,  so  long  as  care  is 
taken  to  insure  that  discontinuities  in  the  material  response  are  not  introduced  at 
the  juncture  between  states  of  reversible  and  irreversible  deformation. 

In  general,  data  for  relatively  soft  rocks  such  as  tuff  display  a nonlinear 
response  even  in  the  elastic  stress  regime.  In  addition,  static  and  dynamic  measure- 
ments of  elastic  moduli  usually  differ.  Non-linearities  in  the  moduli  are  included 
in  the  numerical  model  when  sufficient  information  exists  to  define  the  material 
behavior;  however,  in  most  cases  constant  moduli  are  sufficient  to  adequately 
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characterize  the  elastic  response  unless  the  test  bed  response  is  of  particular  interest 
in  the  low  stress  regime. 

Our  numerical  models  have  generally  utilized  dynamic  values  of  elastic  moduli 
determined  from  seismic  data  for  the  test  site  of  interest.  Ultrasonic  measurements 
conducted  on  small  laboratory  samples  usually  yield  appreciably  higher  velocities 
than  are  obtained  seismically  for  the  material  in  situ,  although  values  of  Poisson's 
ratio  based  on  ultrasonic  velocity  measurements  have  generally  supported  values 
derived  from  seismic  velocity  data. 

2.4.5  Plastic  Response 

Plastic  response  in  the  numerical  model  is  based  on  the  usual  formulation  of 
plasticity  theory  in  which  the  total  strain  rate  tensor  is  separated  into  elastic  and 
plastic  components,  where  the  elastic  components  of  strain  rate  are  assumed  to 
satisfy  Hooke's  law.  Plastic  flow  occurs  only  after  the  stress  state  reaches  a critical 
yield  surface  which  in  general  is  a function  of  both  the  current  state  of  stress  and 
of  the  past  incremental  strain  history.  Thermal  softening  in  response  to  deformation 
heating,  work  hardening,  and  changes  in  effective  stress  due  to  dilatation  or  com- 
paction can  all  potentially  contribute  strain  and  strain  rate  dependent  terms  to  the 
yield  surface. 

Stress  deviators  initially  calculated  on  the  basis  of  Hooke's  law  to  exceed  the 
yield  surface  are  adjusted  to  the  yield  surface  locus  using  a particular  flow  rule. 
The  associative  flow  rule  adjusts  the  stress  deviators  perpendicular  to  the  yield 
surface,  resulting  in  an  increased  mean  stress  (or  equivalently,  material  dilatation) 
if  plastic  flow  occurs  in  a region  where  the  yield  strength  is  increasing  with  confining 
pressure  (such  as  is  the  case  for  a Mohr-Coulomb  surface).  Since  the  mean  normal 
stress  varies  during  the  process  of  adjusting  the  stress  deviators,  iteration  or  a Taylor 
series  approximation  to  the  strength  function  is  required  if  the  associative  flow  rule 
is  invoked.  Alternatively,  adjusting  the  deviators  at  constant  pressure,  the  non- 
associative  flow  rule,  does  not  alter  the  mean  normal  stress,  nor  does  it  introduce 
a plastic  flow  related  dilatation. 

In  addition  to  depending  on  the  shape  of  the  yield  surface,  bulking  produced 
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by  the  associative  flow  rule  depends  strongly  on  the  value  of  Poisson's  ratio,  being 
a maximum  for  v = o and  zero  for  v = 0.5.  Thus  we  generally  avoid  this  formulation 
unless  sufficient  material  response  data  exists  to  support  its  use.  Both  pancake  and 
chimney  calculations  invoked  the  non-associative  flow  rule  to  adjust  the  stress 
deviators  during  plastic  deformation. 

2.4.6  Fracture  Model 

A variety  of  tensile  fracture  models  have  been,  or  currently  are  in  use  in 
(9) 

ground  motion  codes.  They  vary  in  complexity  from  simple  models  to  rather 

sophisticated  time-dependent  theories.  Unfortunately,  the  data  required  to  establish 

values  for  the  constants  in  the  more  complex  theories  are  almost  entirely  lacking 

for  brittle  fracture  in  earthen  media.  Although  some  data  supporting  time-dependent 

fracture  initiation  exists,  they  suggest  that  rate-dependent  fracture  models  are  not 

needed  for  ground  shock  calculations  in  which  time  scales  of  milliseconds  to  seconds 
(9) 

are  important. 

Currently,  the  failure  criterion  in  STAR  is  based  on  the  mean  normal  stress 
reaching  a critical  minimum  value.  In  the  chimney  and  pancake  calculations  discussed 
in  this  report  this  critical  value  was  assumed  to  be  zero.  Whenever  the  mean  normal 
stress  was  calculated  to  fall  below  this  value,  both  the  mean  normal  stress  and  the 
deviators  were  set  to  zero.  If  this  failed  material  subsequently  recompacted  to 
normal  density,  it  was  assumed  to  be  fully  healed  and  thus  able  to  achieve  its  former 
strength.  It  should  be  noted  that  this  failure  model  will  permit  the  material  to 
develop  some  tensile  strength  if  one  or  more  of  the  principal  stress  components  are 
sufficiently  compressive  to  result  in  a net  positive  mean  normal  stress. 

2.5  OVERLAY  TECHNIQUES 

In  order  to  obtain  a more  accurate  solution  at  early-times  and  to  avoid 
problems  with  early-time  distortion  due  to  cavity  growth,  it  has  been  found  convenient 
to  overlay  the  solution  of  one-dimensional  shock  hydrodynamic  problem  into  the 
two-dimensional  STAR  grid.  The  technique  for  accomplishing  this  in  a way  which 
minimizes  the  perturbations  due  to  the  overlay  is  fairly  complicated;  however,  an 
outline  of  the  technique  will  be  given  here. 

Each  cell  in  the  STAR  grid  is  divided  up  into  an  even,  but  otherwise  arbitrary 
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number  of  triangular  subregions  (72  subregions  have  been  used  for  the  present 
calculations).  The  areas  and  centers  of  each  of  these  subregions  are  established  and 
the  interpolated  velocities  for  the  regions  are  computed.  The  following  area-weighted 
variables  for  each  STAR  cell  are  then  summed:  density,  internal  energy,  kinetic 
energy,  momentum,  maximum  density  (crushup  information),  three  of  the  four  stress 
deviators,  and  pressure.  The  fourth  stress  deviator  is  determined  from  the  constraint 
that  the  trace  of  the  deviatoric  stress  tensor  must  vanish.  The  stresses  are  then 
rotated  into  the  proper  coordinate  system  through  the  appropriate  tensor  transfor- 
mation. The  co-areas  are  checked  to  see  if  any  STAR  cell  lies  outside  the  boundaries 
of  the  SIMONE  calculation  and  the  summations  are  adjusted  as  required. 

Density,  mass,  total  energy,  plastic  energy,  minimum  cell  volume  and  the 
three  deviators  are  then  calculated  and  the  equation  of  state  is  entered  with  this 
density  and  energy  to  return  a pressure.  If  the  pressure  returned  is  within  5%  of 
the  value  overlayed  from  the  ID  computation,  it  is  accepted;  otherwise,  the  density 
and  energy  are  established  from  an  iteration  with  the  Hugoniot  until  the  pressure 
returned  by  the  equation  of  state  is  equal  to  the  pressure  from  the  overlay  within 
an  acceptable  bound.  The  velocities  of  the  STAR  grid  intersections  are  then 
interpolated  from  their  location  with  respect  to  the  SIMONE  grid. 

Unfortunately,  this  version  of  the  overlay  technique  was  established  following 
the  computation  of  the  strong  chimney  material  cases.  The  iteration  utilizing  the 
Hugoniot  was  not  used  for  those  cases,  and  as  a result  certain  perturbations  to  this 
solution  resulted.  These  perturbations  affect  the  region  immediately  adjacent  to  the 
cavity  and  will  be  considered  in  the  analysis  in  Section  3. 
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3.  PANCAKE  SIMULATION  OF  NEARBY  CHIMNEY 

3.1  CONFIGURATIONS  OF  PANCAKE  COMPUTATIONS 

3.1.1  Geometry 

The  region  modeled  in  the  axisymmetric  computations  reported  here  included 
locations  extending  to  1200  meters  from  the  working  point  in  the  positive  and 
negative  axial  directions,  as  well  as  in  the  radial  direction.  The  sketch  in  Figure  1 
shows  this  region  and  the  coordinate  system  used  for  the  computations.  The  entire 
mesh  consisted  of  7,938  cells.  Although  this  number  of  cells  would  suggest  relatively 
crude  zoning,  it  is  adequate  to  resolve  the  major  features  of  interest  regarding  the 
interaction  of  the  ground  shock  with  the  chimney  region.  Advantage  has  been  taken 
of  so-called  "sigma  zoning"  wherein  zone  sizes  are  increased  as  a function  of  distance 
away  from  the  origin. 

The  region  included  in  the  computation  has  been  chosen  large  enough  such 
that  waves  do  not  reach  the  boundaries  during  the  time  of  interest  of  the  computation. 
This  completely  obviates  the  necessity  for  sophisticated  boundary  treatments  which 
sometimes  can  impact  solutions  adversely. 

A limited  region  of  the  calculational  grid  which  was  used  for  analysis  is  shown 
in  Figure  2.  The  effects  of  an  explosion-produced  chimney  was  modeled  in  these 
two-dimensional  computations  as  a disk  or  pancake  where  the  area  of  the  pancake 
which  intersects  direct  rays  from  the  working  point  of  the  explosion  was  taken  to 
be  equal  to  the  area  subtended  by  the  expected  chimney  profile  for  a chimney 
centered  at  a range  of  152  meters  (500  feet)  from  the  working  point.  The  thickness 
of  the  pancake  was  selected  to  give  a volume  equal  to  actual  average  chimney 
volumes. 


Two  chimney  positions  were  studied:  a "far"  location  with  the  chimney  centered 
at  152  meters  (500  feet)  from  the  working  point  and  a "near"  chimney  at  a range 
of  122  meters  (400  feet).  The  locations  of  the  "near"  and  "far"  chimneys  are  shown 
in  Figure  2.  The  doubly  cross-hatched  region  indicates  an  overlap  of  pancake  volumes 
from  one  case  to  the  next.  Notice  that  for  the  disk  centered  at  122  meters,  its 
interface  with  the  tuff  surrounding  the  cavity  is  only  95  meters  from  the  working 
point.  One  set  of  computations  was  done  for  a chimney  material  which  was  relatively 
strong  and  another  set  of  computations  was  performed  for  a chimney  material  which 


supported  no  shear  stress  (weak  case).  Prior  to  the  performance  of  the  weak  chimney 
material  set  of  computations,  a technique  to  distort  the  initial  grid  to  reduce  errors 
in  the  overlaying  procedure  (discussed  in  Section  3.1.2)  was  developed.  This  technique 
led  to  the  grid  shown  in  Figure  3 which  replaced  the  equivalent  region  of  the  grid 
shown  in  Figure  2.  Details  of  the  material  properties  of  the  chimneys,  as  well  as 
the  free-field  tuff,  are  discussed  in  Section  3.2  below. 

The  limited  grid  of  Figure  2 will  be  used  as  the  region  in  which  certain 
contour  and  isometric  plots  are  drawn.  Also,  the  darkened  cells  at  about  80  meters 
from  the  working  point  indicate  the  regions  where  time  history  data  has  been  saved. 

3.1.2  Initial  Conditions 

The  two-dimensional  STAR  calculations  performed  were  initialized  using  the 
results  of  a one-dimensional  calculation  performed  with  the  ID  SIMONE  code.  The 
utilization  of  the  ID  code  to  compute  the  early-time  shock  hydrodynamics  is  at  once 
more  accurate  and  less  expensive  than  performing  the  early  dynamics  in  the  2D 
code.  The  large  distortions  which  occur  adjacent  to  the  expanding  cavity  are  more 
easily  accommodated  in  the  ID  code.  Extensive  rezoning  which  would  be  required 
in  the  2D  computation  is  known  to  lead  to  some  diffusion  which  may  impact  the 
solution. 

The  initial  conditions  for  all  four  studies  were  taken  from  a SIMONE  calculation 
of  a hypothetical  10  kT  explosion  in  the  same  tuff  which  was  used  as  the  medium 
for  the  two-dimensional  calculations  (see  Section  3.2  for  properties).  Since  until  the 
wave  reaches  the  inhomogeneity  caused  by  the  inclusion  of  the  pancakes,  the  problem 
is  that  of  a spherical  burst,  it  is,  in  principle,  possible  to  allow  the  initialization  of 
the  problems  for  the  "far"  pancakes  to  occur  at  a later  time  than  for  the  "near" 
pancakes.  A comparison  of  the  early  behavior  for  problems  initialized  at  different 
times  indicated  that  some  differences  did  exist  in  the  solution  in  the  region  which 
was  unaffected  by  either  chimney.  As  a result,  all  four  computations  were  initialized 
at  the  time  of  33.7  milliseconds  in  order  to  improve  comparisons  among  the 
calculations. 

The  overlays  for  the  two  sets  of  computations  were  performed  using  different 
techniques.  The  two  major  differences  in  these  techniques  concerned  one  aspect  of 
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the  overlay  procedure  itself  and  the  utilization  of  a deformed  grid  in  the  two- 
dimensional  problems.  The  original  overlay  procedure  which  was  used  for  the  strong 
pancake  cases  conserved  momentum  and  energy  very  precisely,  but  resulted  in  some 
zones  for  which  the  overlaid  variables  were  not  consistent  with  the  equation-of-state 
of  the  material.  The  weak  pancake  cases  utilized  the  technique  which  is  described 
in  Section  2.5.  The  utilization  of  the  deformed  grid  for  the  weak  pancake  problems 
allowed  us  to  avoid  certain  problems  in  determining  the  co-areas  for  the  "staircase" 
grid  in  the  cavity  regions  shown  in  Figure  2. 

A comparison  of  the  velocity  profiles  in  the  2D  grids  for  the  two  sets  of 
problems  at  a time  of  33.73  ms  is  shown  in  Figure  4.  No  substantial  differences 
are  seen  to  exist  between  the  two  sets  of  problems.  The  axial,  mean,  and  hoop 
stress  profiles  for  the  two  studies  are  shown  in  Figures  5a  and  b.  For  the  stresses, 
some  differences  can  be  observed  near  the  cavity  and  at  the  shock  front.  As 
indicated  earlier,  the  overlay  procedure  which  was  used  reduces  the  perturbations  in 
the  solution  from  incompatabilities  of  overlaid  conditions  with  the  material  equa- 
tion-of-state. 

While  a number  of  test  calculations  have  been  done  to  develop  the  current 
grid  distortion  and  overlay  procedures,  the  extensive  studies  required  to  sort  out 
differences  in  the  computations  due  to  these  modifications  have  not  been  done. 
However,  some  conclusions  can  be  drawn  from  the  solutions  at  hand  and  these  will 
be  discussed  in  Section  3.3. 

3.2  MATERIAL  PROPERTIES  FOR  THE  PANCAKE  SIMULATIONS 
3.2.1  Undisturbed  Tuff 

Three  different  material  regions  reside  in  the  pancake  calculations:  the  cavity, 
the  disturbed  chimney  material,  and  the  surrounding  undisturbed  tuff.  As  discussed 
in  Section  2.4.2,  the  cavity  gas  was  described  by  a simple  gamma-law  gas  treatment. 
Both  the  disturbed  chimney  material  and  the  surrounding  tuff  were  initially  assumed 
to  be  loaded  hydrostatically  to  a uniform  pressure  of  6.83  MPa,  a value  typical  of 
the  lithostatic  overburden  pressure  in  tuff  at  a depth  of  about  400  meters.  A single 
material  response  model  was  employed  for  each  separate  material.  Thus,  layering 
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Figure  4.  Velocity  profile  overlaid  into  the  2D  grid  at  33.73  ms  for  the  pancake  problems 
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Figure  5.  Stress  profiles  overlaid  into  pancake  problems. 
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effects  in  either  the  chimney  or  surrounding  tuff  were  not  considered.  Material 
properties  of  the  undisturbed  tuff  were  selected  to  be  reasonably  representative  of 
the  response  of  typical  NTS  tuff,  although  there  are  deviations  as  will  be  discussed. 
Development  of  the  pancake  material  models  is  described  in  Section  3.2.2. 

Material  models  used  to  described  the  response  of  the  undisturbed  tuff  were 
developed  from  laboratory  test  data  obtained  for  samples  taken  from  various  locations 
in  the  Mighty  Epic  test  bed/10^  Elastic  properties,  listed  in  Table  1,  are  based  on 
USGS  seismic  log  data  obtained  in  exploratory  hole  U12n.l0  UG#7^^in  addition  to 
the  laboratory  response  data  contained  in  Reference  12.  The  resulting  longitudinal 
velocity,  V , of  2286  m/s  is  not  unusually  low,  although  it  is  about  270  m/s  less 
than  the  average  value  derived  from  seismic  data  for  other  test  sites  in  competent 

3 

tuff.  Density  of  the  undisturbed  tuff  prior  to  shock  compression,  1.920  Mg/m  , is 
typical  of  values  measured  for  NTS  tuff. 

Response  of  the  fully  compacted  material  was  represented  in  the  code  by  a 
tabular  equation  of  state  generated  using  the  techniques  discussed  in  Section  2.4.2, 
assuming  a water  content  of  17%  by  weight,  with  dry  rock  comprising  the  balance 
of  the  material.  The  equation  of  state  for  dry  rock  was  calculated  using  a 
Mie-Gruneisen  representation 

P ■ G0Poe  + Ph<v)[1  - -T2-  (,o  - v>]  <3> 

3 

where  Gq  = 0.33  is  the  Gruneisen  ratio,  vQ  = l/po  = 0.4219  m /Mg  is  the  ambient 
specific  volume,  and  pH  is  the  Hugoniot  shock  pressure  derived  from  a quadratic  fit 
of  shock  velocity  to  measured  particle  velocity  data  given  by 


u = A + B u + C u 2 
s P P 

for  constants  A,  B and  C equal  to  350.25  m/s,  0.70476  and  1.0055 
respectively. 


(4) 

x 10-4  s/m, 
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Table  1 

Low  Pressure  Properties  Materials  in  Pancake  Simulations 


Material 

Property 

Undisturbed 

Tuff 

Strong 

Pancake 

Weak 

Pancake 

3 

Density,  pQ  (Mg/m  ) 

1.920 

1.470 

1.470 

Young's  Modulus,  EQ(GPa) 

6.78 

1.54 

0 

Bulk  Modulus,  Kq  (GPa) 

6.64 

0.214 

0.214 

Rigidity  Modulus,  Gq  (GPa) 

2.55 

2.55 

0 

Compressional  Wave  Speed,  Cp(m/s) 

2286 

1568 

382 

Shear  Wave  Speed,  Cg  (m/s) 

1152 

1317 

0 

Water  Content  by  Weight,  Mw(%) 

17 

17 

17 
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Irreversible  void  collapse  at  pressures  in  excess  of  the  elastic  limit  was  modeled 
using  the  crush  response  illustrated  in  Figure  6,  where  the  void  collapse  for  the 
pancake  rubble  is  also  shown  for  comparison.  The  initial  2%  void  content  for  the 
undisturbed  tuff  is  somewhat  higher  than  average  values  measured  for  recent  tests 
(other  than  Mighty  Epic)  which  have  ranged  from  about  1.2  to  1.6%.  However,  a 
value  of  2%  air  voids  is  well  within  the  normal  spread  of  individual  sample  values 
which  range  from  less  than  1%  to  over  5%,  and  2%  is  far  less  than  average  values 
obtained  for  some  older  tests  conducted  in  unsaturated  media.  The  procedure 
described  in  Section  2.4.3  was  used  to  specify  release  paths  from  the  crush  curves 
of  Figure  6.  Constants  used  to  specify  the  shape  of  the  crush  curve  are  summarized 
in  Table  2. 

The  pressure  and  energy  dependent  yield  function,  Y,  employed  in  the  calcu- 
lations is 

Y = Yn  + Y E ( 2 • E I / 1 - I \ ; 

m \ ^m/  \ em  / 

Y ' (¥o  + VJ('  - tj  ; 

Y = 0 

where  Yo  is  the  unconfined  yield,  Ym  is  the  maximum  increase  in  yield  strength, 
P is  the  pressure  at  which  maximum  strength  is  attained,  and  e_  is  the  melt 
energy  of  the  material,  taken  to  be  1.91  x 10  J.  Values  of  Pm,  Yq  and  Ym  are 
summarized  in  Table  2 for  the  various  pancake  simulation  materials,  and  the  resulting 
yield  surface  is  illustrated  in  Figure  7. 

3.2.2  Chimney  Materials 

At  the  time  these  studies  were  initiated  very  little  data  other  than  qualitative 
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Chimney  Material 


Figure  6 
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Porous  crushup  behavior  for  undisturbed  tuff  and  chimney  material 
in  pancake  computaions. 


Table  2 

Crush  Curve  and  Yield  Surface  Constants  in  Pancake  Simulations 


Material 

Property 

Undisturbed 

Tuff 

Strong 

Pancake 

Weak 

Pancake 

Ambient  Density,  p0(Mg/m3) 

1.92 

1.47 

1.47 

Ambient  Solid  Density,  ps(Mg/m3) 

1.96 

1.96 

1.96 

Air  Void  Content,  V (%) 

fi 

2.0 

25.0 

25.0 

Elastic  Pressure,  p (MPa) 

0 

10 

10 

10 

Crush  Pressure,  pc  (MPa) 

400 

400 

400 

Unconfined  Yield,  Yq  (MPa) 

18.0 

18.0 

0 

Yield  Strength  Increase,  Ym  (MPa) 

52.0 

52.0 

0 

Maximum  Yield  Increase 

Pressure,  pm  (MPa) 

200 

200 

0 

observations  of  chimney  rubble  were  available  to  characterize  the  response  of  such 
materials  to  ground  shock  loads.  This  difficulty  led  to  the  two  sets  of  bounding 
calculations  which  assume  rather  extreme  limits  for  the  strength  of  the  chimney 
material. 

In  the  pancake  simulations  the  chimney  rubble  was  assumed  to  have  the  same 
composition  as  surrounding  tuff.  Thus,  the  tabular  equation  of  state  described  in 
Section  3.2.1  was  also  utilized  to  specify  the  response  of  the  chimney  material  in 
the  fully  crushed  state.  Of  course,  the  initial  air  void  content  and  crush  curve 
response  of  the  pancake  material  were  much  different  from  that  of  the  surrounding 
tuff  (see  Table  2 and  Figure  6).  The  initial  air  void  content  of  the  chimney  material 
was  derived  by  assuming  that  the  cavity  volume  produced  by  a nominal  10  kt  explosion 
would  subsequently  be  distributed  uniformly  throughout  the  chimney  rubble  following 
chimney  collapse.  The  chimney  volume  used  was  based  on  limited  drill-back  data 
which  suggest  an  average  chimney  radius  of  1.1  times  the  final  cavity  radius,  Rc, 
and  total  height  of  5.13  times  Rc.  A cylindrical  chimney  shape  with  hemispherical 
ends  was  assumed  since  little  data  were  available  regarding  shape  of  the  chimney 
above  tunnel  level. 

The  crush  curve  shown  in  Figure  6 for  pancake  material  is  relatively  uncertain. 

(13) 

It  was  derived  from  uniaxial  test  data  obtained  in  exploratory  hole  UE12nlO#9. 

The  crush  response  determined  for  the  sample  having  the  highest  air  void  content, 
10.3%,  was  scaled  linearly  in  specific  volume  to  reflect  the  25%  void  content  required 
by  distributing  the  chimney  volume  throughout  the  chimney  material  (see  Figure  6). 
The  initial,  elastic  bulk  modulus  was  chosen  to  be  compatible  with  the  zero  pressure 
slope  of  the  crush  curve;  whereas,  in  the  strong  chimney  case  the  rigidity  modulus 
was  assumed  to  be  the  same  as  the  undisturbed  tuff,  in  keeping  with  the  assumption 
that  the  failure  response  was  unaltered.  For  the  zero  strength  pancake  simulation 
the  rigidity  modulus  was  set  to  zero  since  deviatoric  stresses  could  not  be  supported. 

3.3  RESULTS 

3.3.1  Strong  Chimney  Material 

The  ground  shock  from  the  detonation  propagates  outward  toward  the  pancake 


regions  until  rarefactions  are  produced  at  about  42  ms  from  the  near  pancake; 
however,  the  wave  does  not  reach  the  far  pancake  until  55  ms.  By  60  ms  the 
rarefaction  from  the  near  pancake  has  propagated  back  to  a range  of  50  meters 
from  the  detonation  point  producing  very  little  effect  at  that  range,  but  causing 
about  a 50  percent  reduction  in  axial  stress  near  the  chimney  boundary  itself  (Figure 
8).  At  this  time,  the  effects  due  to  the  far  pancake  are  nil. 

The  rarefaction  from  the  near  pancake  continues  to  impact  the  axial  stress 
as  time  goes  on,  but  at  50  meters  at  a time  of  80  ms  the  difference  between  the 
free-field  and  the  near  pancake  axial  stress  is  down  to  about  20  percent  (Figure  9). 
By  80  ms  the  effect  of  the  far  pancake  has  propagated  back  to  a range  of  100 
meters,  but  its  effect  is  seen  to  be  less  than  20  percent  at  about  120  meters  which 
is  immediately  in  front  of  the  far  chimney. 

The  displacements  of  points  initially  80  meters  from  the  working  point  are 
shown  in  Figure  10.  The  displacement  response  for  the  far  pancake  problem  is 
essentially  the  same  as  that  for  the  free-field.  It  is  seen  that  the  maximum 
difference  in  displacement  for  the  near  pancake  occurs  at  about  125  ms  and  is  only 
0.22  meters.  This  displacement  difference  endures  until  250  ms  when  the  solution 
was  ended.  At  this  time  velocities  were  approximately  zero  throughout  the  region 
of  interest  with  only  a small  elastic  wave  propagated  to  the  far  field. 

As  indicated  above,  axial  stress  has  been  considered  in  order  to  evaluate  the 
propensity  for  dynamic  cracking  to  occur  between  the  chimney  and  the  detonation 
point  as  a result  of  the  rarefaction  wave  emanating  from  the  chimney.  High  values 
of  hoop  and  radial  stresses  are  believed  to  be  responsible  for  preventing  the  late-time 
hydrofracture  of  the  site  medium  by  cavity  gases.  A hoop  stress  time  history  for 
the  field  point  at  80  meters  is  shown  in  Figure  11.  The  far  pancake  history  is  very 
close  to  that  of  the  free-field,  whereas  at  50  ms,  the  close  pancake  solution  diverges 
from  the  free-field  response,  exhibiting  tension  at  around  60  ms.  Recall  that  the 
axial  stress  was  also  depressed  with  respect  to  the  free-field  values  at  this  time 
and  at  a range  of  80  meters  (See  Figure  8).  By  90  ms,  however,  the  near  pancake 
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acement  of  a point  initially  80  m from  the  WP. 


stress  at 


tangential  stress  is  very  close  to  that  of  the  free-field.  During  the  next  reverberation 
it  becomes  depressed  oice  again,  however,  the  rebound  which  commences  at  about 
130  ms  d-ives  the  tangential  stresses  for  all  cases  upward  to  the  rebound  values 
before  the  residual  stresses  establish  themselves.  The  residual  tangential  stresses 
for  the  two  pancake  problems  are  each  within  about  5 percent  of  that  for  the 
free-field.  Considering  that  the  edge  of  the  chimney  is  only  about  15  meters  from 
the  field  point  being  considered  it  appears  that  the  inclusion  of  the  pancake  for 
both  cases  has  produced  relatively  small  effects  in  the  late-time  stress  fields. 

The  residual  hoop  stress  about  the  cavity  at  a time  of  250  ms  is  shown  in 
Figure  12  for  the  near  pancake  computation.  The  predominant  influence  is  seen  to 
be  confined  to  a region  relatively  close  to  the  pancake  itself  and  outside  the  region 
of  high  residual  stress  around  the  cavity.  A much  better  feeling  for  the  influence 
of  the  chimney  on  the  stress  field  can  be  obtained  from  Figure  13,  in  which  the 
contour  plot  of  the  close  chimney  hoop  stress  field  (dashed  lines)  has  been  superimposed 
on  the  results  obtained  for  the  farther  chimney  location  (solid  contour  lines).  Although 
the  stress  field  is  modified  slightly  near  the  cavity  in  the  two  cases,  the  region  in 
which  the  stress  field  is  significantly  affected  is  seen  to  lie  well  beyond  the  peak 
in  the  residual  stress  field.  It  is  interesting  to  note  that  the  closer  chimney  location 
appears  to  increase  the  peak  hoop  stress  very  slightly  rather  than  to  degrade  it  in 
this  region  as  one  might  expect.  Of  course,  in  the  vicinity  of  the  chimney  the 
residual  hoop  stress  is  reduced  appreciably  by  the  stress  relief  in  the  weak  rubble 
material. 

The  tangential  stresses  for  both  cases  are  shown  for  the  ray  connecting  the 
detonation  point  with  the  center  of  the  pancake  in  Figure  14.  The  influence  of 
even  the  close  chimney  is  seen  to  be  negligible  at  ranges  less  than  70  meters.  Even 
out  to  the  edge  of  the  near  chimney,  the  differences  are  not  substantial,  being  a 
maximum  of  about  15  percent. 

The  residual  stress  fields  for  the  far  pancake  problem  are  shown  in  Figure 
15.  In  this  important  region,  adjacent  to  the  cavity,  the  shape  of  the  stress  fields 
do  not  show  any  unusual  characteristics  when  compared  with  results  from  other 
studies. 
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Figure  13 


Comparison  of  tangential  stress  near  the 
cavity  for  the  two  pancake  problems  at 
250  ms.  Contour  interval  is  5 MPa; 
solid  contours  are  for  Far  chimney,  dashed 
contours  denote  results  for  close  chimney 
location.  Plus  symbol  marks  every  fifth 
contour. 
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Figure  15.  Tangential  and  axial  stress  at  250  ms  calculated  along  the  axis 
of  symmetry  for  the  pancake  located  152.5  m (500  ft)  from  the  WP 


3.3.2  Weak  Chimney  Material 

The  results  from  the  strong  chimney  compuations  discussed  above  indicated 
that  the  far  chimney  had  very  little  influence  on  the  stress  field  around  the  cavity 
at  any  time.  The  near  chimney  (centered  at  122  meters)  produced  a rarefaction 
wave  which  was  strong  enough  to  result  in  a modest  and  temporary  reduction  of  the 
stresses  in  the  vicinity  of  the  cavity  at  a time  of  about  80  ms.  Neither  chimney 
location  produced  more  than  a slight  change  in  the  residual  stress  field  around  the 
cavity. 


There  was  concern  that  the  material  model  for  the  chimney  rubble  used  in 
these  first  results  might  have  been  too  strong;  thus,  underestimating  the  possible 
effect  on  the  shock  wave  interaction  and  residual  stress  field.  In  order  to  establish 
some  bound  on  the  effects  of  this  uncertainty,  the  strength  of  the  rubble  was  reduced 
to  zero  in  the  calculations  here.  This  is  clearly  overly  conservative,  particularly  at 
shot  level  in  the  chimney  where  drill-back  indicates  that  the  strength  of  the  material 
within  the  chimney  is  comparable  to  that  of  the  surrounding  tuff. 

The  effect  of  the  rarefaction  from  the  chimney  region  on  the  axial  stress 

profiles  at  times  of  60,  80  and  100  ms  is  shown  in  Figures  16  through  18,  respectively. 

The  rarefaction  has  reached  a range  of  approximately  45  meters  at  60  ms  (Figure  16) 

in  the  near  chimney  calculation.  The  dramatic  reduction  of  axial  stress  adjacent 

to  the  weak  pancake  is  already  apparent  at  this  time.  An  effective  spall  into  the 

chimney  reduces  the  axial  stress  immediately  in  front  of  the  chimney  to  zero.  This 

type  of  separation  is  also  apparent  in  stemming  computations  at  interfaces  between 

(6) 

a rock  and  weak  stemming  material. 

At  60  ms  the  wave  is  just  reaching  the  close  edge  of  the  far  pancake;  hence, 
the  axial  stresses  in  the  problem  are  not  yet  perturbed  by  its  presence.  The  far 
chimney  axial  stress  profile  is  identical  to  that  of  the  free-field  (taken  as  the  ray 
away  from  the  near  chimney.) 

The  hoop  stress  profile  at  60  ms  (Figure  19)  shows  that  the  rarefaction  leads 
to  tension  in  this  component  from  the  pancake  edge  at  95  meters  back  to  65  meters 
from  the  detonation  point.  The  radial  stress  component  is  also  tensile  in  this  region, 
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(Note  suppressed  zero  on  the  abscissa. 


calculations.  (Note  suppressed  zero  on  the  abscissa. 


but  is  not  identical  to  the  hoop  stress  due  to  the  asymmetry  produced  by  the  disk. 
At  times  prior  to  the  rarefaction,  the  hoop  and  radial  stresses  are  essentially  identical. 
The  far  pancake  profile  is  unperturbed  at  this  time. 

By  80  ms  a large  region  has  been  fractured  between  75  meters,  and  the 
pancake  edge  and  the  peak  axial  stress  is  reduced  by  nearly  a factor  of  two  by  the 
presence  of  the  near  chimney  (Figure  17).  Also,  the  combination  of  the  rarefaction 
from  the  cavity  and  that  from  the  pancake  have  caused  the  unloading  of  the  tuff 
adjacent  to  the  cavity.  The  far  chimney  has  caused  a rarefaction  which  has  propagated 
back  to  60  meters.  This  wave  has  not  yet  influenced  the  peak  axial  stress  near  the 
cavity,  but  has  created  a 35  percent  drop  in  stress  at  about  100  meters. 

The  region  in  the  near  pancake  problem  which  was  fractured  at  80  ms  has 
started  to  reload  by  100  ms  (Figure  18).  The  rock  located  between  69  and  79  meters 
is  fractured  (in  the  sense  that  this  material  is  below  normal  density  and  mean  stress 
is  zero)  while  rock  closer  to  the  chimney  edge  is  alternately  in  tension  and 
compression.  While  it  is  believed  that  reloading  is  probable  back  to  the  overburden 
level  it  should  be  borne  in  mind  that  the  numerical  treatment  of  fracture  may  play 
a role  in  establishing  the  stress  state,  especially  in  the  immediate  vicinity  of  a 
fractured  region. 

The  tangential  stress  profile  for  the  near  pancake  at  100  ms  shows  the 
fractured  region  between  69  and  79  meters  (Figure  20).  The  peak  hoop  stress  is 
about  the  same  as  the  free-field  value  as  was  the  case  for  the  axial  stress,  however 
the  region  over  which  high  stress  acts  is  about  one  half  of  that  for  the  free-field 
case.  The  reloading  of  the  hoop  stress  near  the  chimney  to  levels  above  overburden 
is  credible  as  an  arching  mechanism  around  the  weak  chimney.  Similar  effects  have 
been  noted  in  computations  of  shock  interactions  with  weak  sand  columns/®) 

For  the  far  chimney  at  100  ms  no  region  of  fracture  exists.  The  axial  stresses 
for  the  10  meters  in  front  of  the  pancake  are  seen  to  be  approaching  zero  (Figure  18). 
The  rarefaction  has  influenced  the  axial  stress  back  to  the  cavity  by  this  time 
reducing  the  peak  stress  by  only  about  20  percent.  The  hoop  stresses  of  this  time 
are  tensile  near  the  pancake  and  the  peak  is  down  by  approximately  40  percent 


Away  from  Near  Chimney 


(Figure  20).  Arching  around  the  far  pancake  is  also  evident  in  the  figure. 

The  influence  of  the  chimney  at  a point  initially  located  between  the  chimney 
and  the  WP  at  a range  of  80  meters  from  the  WP  is  illustrated  in  Figures  21  and 
22.  The  more  distant  chimney  produces  a relatively  small  change  in  displacement 
( ~20  percent)  and  peak  hoop  stress.  The  close  chimney  location  increases  the  final 
displacement  by  nearly  a factor  of  three  and  the  hoop  stress  is  reduced  to  zero. 

The  large  displacements  are  consistent  with  the  unloading  of  this  region  in  the  near 
pancake  problem.  For  the  distant  chimney  location,  neither  the  dynamic,  nor  residual 
stresses  are  appreciably  influenced  (Figure  22);  however,  the  perturbation  at  this  80 
meter  range  which  results  from  the  presence  of  the  near  chimney  is  evident.  In 
this  case,  the  material  fractures  shortly  after  the  initial  shock  passage  as  the  tuff 
at  the  chimney  edge  spalls  into  the  weak  chimney  region. 

Residual  stresses  are  established  after  rebound  when  only  small  amplitude 
elastic  oscillations  continue  out  far  from  the  source.  The  axial  stress  profiles  (Figure 
23)  show  the  effects  of  the  zero  strength  pancakes  dramatically.  The  near  pancake 

* 8 

has  reduced  the  peak  axial  residual  stress  by  30  percent  and  narrowed  the  region 
near  the  cavity  over  which  this  stress  is  compressive  by  more  than  a factor  of  three. 

The  far  chimney  produces  only  a 15  percent  reduction  of  peak  stress  with  compressive 
stresses  enduring  out  to  87  meters. 

Both  chimneys  leave  fractured  regions  on  axis.  These  regions  are  12  meters 
long  for  the  near  and  15  meters  for  the  far  pancake.  The  near  pancake  fracture 
zone  starts  as  close  as  69  meters  from  the  detonation  point  while  that  for  the  far 
pancake  starts  at  87  meters.  The  hoop  stresses  (Figure  24)  do  not  show  a reduction 
of  the  peaks  as  dramatically  as  that  for  the  axial  stress  components,  but  the  steep 
gradients  to  the  fractured  regions  are  similar.  One  does  not  expect  as  large  a 
reduction  in  the  hoop  stresses  as  in  the  axial  stresses  near  the  cavity  since  the 
rarefaction  from  the  finite  diameter  pancakes  suffers  geometric  attenuation  as  it 
propagates  back  toward  the  cavity. 

3.3.3  Strong  Versus  Weak  Cases 

As  indicated  in  Section  3.1.2,  the  two  sets  of  problems  were  initialized 
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TIME  (ms) 

Figure  21.  Comparison  of  displacement  of  points  initially  80  meters  from  WP 


TIME  (ms) 

Figure  22.  Tangential  (hoop)  stress  at  points  initially  80  meters  from  the  WP 


Figure  24.  Comparison  of  tangential  stress  profiles  at  241  ms  in  the  weak  chimney 
calculations.  (Note  suppressed  zero  in  the  abscissa.) 


differently.  The  differences  were  in  the  utilization  of  a deformed  grid  (Figure  3) 
in  the  vicinity  of  the  cavity  and  an  overlay  technique  which  produced  self-consistent 
field  variables  for  the  weak  cases.  The  results  indicate  that  at  times  before  about 
80  ms  it  is  not  possible  to  directly  compare  the  two  sets  of  problems.  For  example, 
the  axial  stress  profiles  at  60  ms  for  the  differing  material  cases  (Figures  8 and  16) 
show  a region  from  45  to  65  meters  for  which  the  axial  stress  in  the  strong  chimney 
material  cases  is  higher  than  that  for  the  weak  cases.  Considering  just  the  far 
pancake  problems  which  have  not  been  influenced  by  the  rarefaction  from  the  pancake 
at  this  early  time,  it  is  clear  that  this  difference  must  be  due  to  the  initialization 
of  the  problems. 

While  it's  difficult  to  sort  out  the  precise  cause  for  the  difference,  the  likely 
source  of  this  effect  is  due  to  the  difference  between  the  staircase  cavity  boundary 
utilized  in  the  strong  cases  from  the  more  spherical  cavity  boundary  afforded  by 
the  utilization  of  the  deformed  grid  shown  in  Figure  3.  Because  of  the  flattened 
region  of  cavity  along  the  axis  of  the  problem  where  most  of  our  comparisons  have 
been  made,  a plane  wave  is  driven  into  the  2D  grid  when  in  fact  it  should  be  a 
spherical  wave.  The  difference  in  geometric  attenuation  then,  seems  to  be  responsible 
for  this  effect.  Comparison  of  the  axial  stress  profiles  at  80  ms  (Figures  9 and  17) 
shows  that  the  difference  between  the  peaks  located  at  about  50  meters  from  the 
detonation  point  is  only  about  10  percent  at  this  time.  A comparison  with  one- 
dimensional computations  in  similar  materials  which  have  been  run  out  to  later  times 
shows  that  the  stress  profile  for  the  weak  pancake  material  problems  is  the  more 
correct  one  at  60  ms. 

The  most  dramatic  difference  between  the  two  sets  of  problems  is  clearly 
the  generation  of  significantly  reduced  stresses  and  regions  of  fracture  by  the  weak 
pancake  material. 

The  residual  stress  fields  which  were  established  for  the  two  cases  are 
remarkably  similar  considering  the  fact  that  the  weak  pancake  material  problems 
are  considered  to  be  extreme  bound.  The  residual  hoop  stresses  for  all  four 
computations  are  very  similar  out  to  the  region  where  the  weak  pancake  material 
oases  drop  off  into  the  zones  of  fracture  (Figure  24).  A similar  statement  can  be 


made  for  the  axial  residual  stress  which,  of  course,  unloads  to  a greater  extent  due 
to  the  geometry  of  the  problem. 

The  bounding  computations  including  the  zero  strength  material  in  the  pancakes 

were  performed  because  no  information  was  available  regarding  the  properties  of 
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chimney  material.  Data  obtained  more  recently  suggest  that  the  reduction  in 
the  strength  of  chimney  material  is  only  about  25  percent  at  confining  pressures 
above  1 kbar  along  with  a reduction  of  the  cohesion  to  zero.  This  reduction,  which 
under  confinement  is  within  the  variation  of  the  strength  of  the  rock  from  place  to 
place  in  situ,  is  more  compatible  with  the  strong  pancake  material  problems  than 
the  zero  strength  problems. 

Additional  effects  which  lead  one  to  believe  that  all  of  these  computations 
are  conservative  are  the  placement  of  the  pancake  with  respect  to  the  axis  and  the 
distribution  of  the  cavity  void  volume  uniformly  throughout  the  pancake.  The  former 
effect  is  conservative  because  in  an  actual  configuration,  the  rarefaction  would  not 
be  focused  back  along  the  axis  since  the  chimney  volume  is  distributed  upwards  and 
at  increasing  ranges  from  the  detonation  point.  This  means  that  considerably  more 
geometric  attenuation  can  be  expected  regarding  the  propagation  of  the  rarefaction 
back  towards  the  vicinity  of  the  cavity.  The  latter  effect  has  the  result  of  increasing 
the  strength  of  the  rarefaction  close  to  the  cavity.  This  is  unrealistic  since  the 
region  of  the  chimney  closest  to  the  cavity  is  known  to  be  recompacted  to  nearly 
normal  densities.  A third  difference  between  the  real  world  and  the  computations 
reported  here  is  that  the  pancake  configuration  includes  a plane  surface  from  which 
a nearly  plane  wave  rarefaction  results.  In  the  real  world,  the  chimney  looks  more 
like  a cylinder  whose  axis  is  perpendicular  to  the  axis  of  symmetry  in  the  computations; 
a configuration  which  would  also  produce  significant  geometric  attenuation  of  the 
rarefaction  as  it  moved  back  towards  the  cavity  region. 


4.  IN-CHIMNEY  DETONATION 


4.1  CONFIGURATIONS  OF  IN-CHIMNEY  COMPUTATIONS 

4.1.1  Geometry 

Two  two-dimensional  computations  were  performed  for  the  in-chimney  detona- 
tions study.  The  biggest  differences  between  these  two  calculations  was  in  the 
presumed  shape  of  the  chimney  region  and  the  material  properties  as  discussed  below 
in  Section  4.2.  One  computation  considered  a chimney  which  was  cylindrical  having 
hemispherical  top  and  bottom.  The  second  calculation  included  a conoidal-shaped 
chimney  as  shown  in  Figure  25.  The  natural  symmetry  axis  for  these  problems  lies 
in  the  vertical  direction.  Other  than  the  inclusion  of  a uniform  and  hydrostatic 
overburden  pressure,  the  effects  of  gravity  were  not  modeled  in  these  computations. 
Also,  for  purposes  of  these  exploratory  computations,  the  incorporation  of  layering 
of  the  material  properties  as  a function  of  elevation  was  not  included.  The 
computational  grid  included  a total  of  7200  cells  for  each  problem  and  covered 
distances  of  1060  meters  and  530  meters  in  axial  and  radial  extent,  respectively. 
One-dimensional  calculations  were  used  to  initialize  both  of  the  two-dimensional 
problems  (see  Subsection  4.1.2  below).  For  comparative  purposes  a one-dimensional 
calculation  was  run  to  late-times. 

The  shape  of  the  cylindrical  chimney  shown  in  Figure  25  was  originally  thought 
to  be  characteristic  of  many  of  the  chimneys  explored  in  tuff.  The  diameter  was 
chosen  based  on  a survey  of  measured  chimney  dimensions/1  ^ These  data  suggest 
that  the  chimney  radius  at  the  working  point  level  averages  5 to  15  percent  larger 
than  the  initial  cavity  radius.  Likewise,  the  height  of  the  chimney  was  specified 
as  5.13  times  the  cavity  radius  based  on  Reference  15.  The  volume  utilized  for  the 
cylindrical  chimney  is  consistent  with  that  used  in  the  pancake  studies  discussed  in 
Section  3. 

The  conoidal  shape  shown  in  Figure  25  was  developed  based  on  new  data 
coming  from  exploration  in  progress  at  about  the  time  these  calculations  were 
performed.  ' By  way  of  comparison,  the  lower  hemispheres  for  both  computations 
are  essentially  the  same  as  is  the  overall  height  from  the  very  top  of  the  chimney 
to  the  bottom  of  the  cavity.  The  major  difference  then  is  the  additional  volume 
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due  to  the  flared  side  walls  of  the  chimney.  This  volume  amounts  to  4.87  x 10  m 

5 3 

of  additional  volume  compared  with  a volume  of  4.57  x 10  m for  the  cylindrical 
chimney.  This  volume  difference  will  have  an  effect  on  the  distribution  of  air  voids, 
as  discussed  in  Section  4.2. 

Regions  of  the  computational  mesh  are  shown  for  each  of  the  problems  in 
Figure  26.  Notice  that  the  deformed  grid  technique  discussed  in  Section  3 has  been 
utilized  for  both  problems  in  the  vicinity  of  the  cavity.  No  attempt  has  been  made 
to  deform  the  grid  to  match  the  shape  of  the  chimney  itself.  Since  the  overlay 
occurs  well  inside  this  boundary,  the  impact  of  the  staircase  grid  at  the  chimney 
boundary  should  have  only  minimal  effect  upon  the  results. 

4.1.2  Initial  Conditions 

Initial  conditions  for  the  in-chimney  study  were  obtained  from  two  one- 
dimensional SIMONE  computations.  While  the  initial  conditions  for  pancake  problems 
(Subsection  3.1.2)  were  all  taken  from  one  SIMONE  calculation,  this  was  not  possible 
for  the  in-chimney  studies  due  to  the  material  differences  in  the  chimney  region. 
As  a result,  the  initialization  process  was  performed  when  the  shock  in  each  of  the 
one-dimensional  calculations  reached  a range  of  29.3  meters.  This  corresponded  to 
a time  of  12.84  ms  for  the  cylindrical  chimney,  but  only  7.69  ms  for  the  conoidal 
chimney  due  to  the  higher  propagation  velocity  in  that  material. 

The  velocity  profiles  overlaid  to  the  two  problems  are  shown  in  Figure  27. 
The  steep  wave  front  which  is  calculated  in  the  fine  zone  one-dimensional  calculation 
must  of  necessity  be  smeared  out  in  the  two-dimensional  calculations  due  to  the 
larger  zone  size  employed.  The  overlaid  stresses  for  the  cylindrical  chimney  problem 
are  shown  in  Figure  28.  The  peaks  of  the  stress  wave  which  occur  at  about  26 
meters  have  been  clipped  by  the  overlay  routine  in  order  to  make  the  intra-cell 
thermodynamic  conditions  consistent.  A similar  situation  exists  for  the  stresses 
overlaid  into  the  conoidal  chimney  problem  (Figure  29). 

4.2  MATERIAL  PROPERTIES  FOR  THE  IN-CHIMNEY  DETONATION 
4.2.1  Undisturbed  Tuff 

The  constitutive  model  described  in  Section  3.2.1  was  used  to  specify  the 
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Figure  26.  Interior  sections  of  computational  grids  for  the  in-chimney 
probl ems . 


RANGE  (m) 

a)  Cylindrical  Chimney,  t = 12.84  ms 


10  14  18  22  26 

RANGE  (m) 


b)  Conoidal  Chimney,  t = 7.69  ms 


Figure  27.  Velocity  profiles  overlaid  into  the  two-dimensional  grid 
for  the  in-chimney  problems. 


Figure  23.  Stress  profiles  overlaid  irto  the  cylindrical  chimney  problem  at  12.84  ms.  Radial  stress 
in  the  outward  direction  from  the  WP. 


the  outward  direction  from  the  WP. 


response  of  the  undisturbed  tuff  in  the  pancake  problems.  Thus  differences  between 
the  ground  motion  calculated  for  the  pancake  and  the  in-chimney  studies  are  not 
due  to  differences  in  the  site  material  properties.  As  discussed  in  Section  3.2.1, 
the  undisturbed  tuff  properties  are  reasonably  representative  of  the  average  response 
measured  for  core  samples  of  Area  12  tuff  at  tunnel  level.  Air  void  content  and 
strength  in  the  calculational  model  are  somewhat  higher  than  average  values  for 
sites  other  than  Mighty  Epic,  while  compressional  and  shear  wave  velocities  are 
slightly  lower.  However,  the  assumed  values  fall  well  within  the  range  of  data 
measured  at  virtually  all  sites  in  competent  tuff. 

4.2.2  Chimney  Materials 

In  the  interim  between  initiating  the  pancake  and  the  in-chimney  calculations, 
additional  material  response  data  became  available  for  chimney  and  "chimney-like" 
materials.  Also,  drill  back  and  chimney  gas  flow  data  were  acquired  during  the 
course  of  the  in-chimney  study  which  assisted  in  locating  the  boundary  of  the 
conoidal-shaped  chimney.  Thus  in  some  cases,  the  material  properties  and  chimney 
volumes  assumed  in  the  in  chimney  computations  differ  in  several  important  respects 
from  the  values  assumed  in  the  pancake  studies  discussed  in  Section  3. 

4.2.2. 1 Cylindrical  Chimney.  As  discussed  in  Section  4.1.1,  chimneys  in  tuff  originally 
were  thought  to  be  cylindrical,  and  this  assumption  was  made  in  computing  the 
chimney  volume  assumed  in  the  pancake  studies  (Section  3.2.2).  This  geometry  was 
also  adopted  for  the  first  in-chimney  calculation.  Distributing  the  cavity  volume 
throughout  the  available  chimney  material  again  led  to  an  initial  air  void  content 
of  25%  in  the  chimney  rubble.  Thus,  for  this  chimney  configuration  the  crush  curve 
utilized  in  the  pancake  simulations  was  adopted  (Figure  6,  Section  3.2.2).  In  accord 
with  our  previous  assumption  that  the  composition  of  the  chimney  rubble  would  be 
the  same  as  the  surrounding  tuff,  the  water  content  of  the  rubble  was  taken  to  be 
17%,  and  the  tabular  equation  of  state  discussed  in  Section  3.2.1  was  used  to  specify 
the  p-v  response  of  fully  crushed  rubble. 

Particular  concern  was  given  to  modeling  the  strength  of  the  chimney  material 
in  view  of  the  results  of  the  pancake  studies  which  amply  demonstrated  the  sensitivity 
of  the  calculated  stress  fields  to  this  property.  Unfortunately,  no  data  were  available 


to  directly  determine  the  strength  of  the  chimney  material  for  events  conducted  in 

Area  12  tuff.  As  a consequence,  we  were  forced  to  model  the  rubble  strength  based 

on  data  obtained  from  rubble  manufactured  by  crushing  samples  of  Tonopah  tuff, 

reconstituting  them  under  a confining  pressure  and  then  measuring  the  strength  of 
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this  "chimney  simulant".  These  data  were  then  compared  to  the  response  of  the 
virgin  material  to  deduce  the  degradation  in  strength  which  would  result  from 
preexisting  fractures.  However,  the  failure  surface  determined  in  this  fashion  could 
not  be  used  directly,  since  Tonopah  tuff  is  a welded  tuff  which  is  far  stronger  than 
tunnel  bed  tuff  in  Area  12.  Instead,  we  utilized  an  average  failure  surface  for  Area 
12  tuff  as  reported  in  Reference  18  scaled  down  as  a function  of  pressure  by  the 
same  ratio  as  was  determined  in  the  Tonopah  tuff  tests.  The  resulting  yield  surface 
is  compared  in  Figure  30  with  those  developed  for  the  surrounding  tuff  and  for  the 
chimney  rubble  in  the  conoidal-shaped  chimney  computation. 

A limited  number  of  tests  have  since  been  conducted  on  virgin,  reconstituted, 
and  pulverized  samples  of  tuff  from  drill  hole  UE12e.l/^  These  data  qualitatively 
support  the  strength  reduction  estimate  developed  above  from  the  Tonophah  data, 
at  least  for  the  7 MPa  confining  stress  level  at  which  data  was  obtained.  More 
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recent  average  tuff  yield  data  in  n-tunnel  generally  have  given  somewhat  greater 
values  for  the  material  strength  than  was  the  case  for  the  average  strength  reported 
in  Reference  18  (see  Figure  30.),  This  suggests  that  the  assumed  rubble  strength 
in  the  cylindrical  chimney  calculations  may  be  too  low,  particularly  for  rubble  at 
and  below  the  tunnel  horizon  where  considerable  compaction  by  the  weight  of  overlying 
material  seems  certain. 

4. 2. 2. 2 Conoidal  Chimney.  The  rather  dramatic  earth  motion  and  residual  stress 
field  obtained  in  the  first  in-chimney  shot  calculation  suggested  that  establishing  a 
satisfactory  residual  stress  field  about  the  detonation  point  might  not  be  possible 
within  the  limits  of  plausible  material  response  models  for  the  chimney  rubble.  To 
address  this  question,  at  least  in  part,  a final  calculation  was  accomplished  which 
incorporated  new  chimney  shape  and  properties  data  which  had  just  become  available. 
Rubble  properties  were  assumed  to  be  the  s^me  as  those  measured  for  the  layer  of 
tuff  which  resides  above  the  Might  Epic  test  bed  tuff,  and  the  numerical  response 
model  was  based  on  seismic  data  and  laboratory  measurement  on  core  samples 
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obtained  from  this  region. 


The  equation  of  state  for  this  material  has  the  form 


(12) 


p$  = Tpe  + Ay$  + gys2  (6) 


O 

where  P,  A and  B are  constants  listed  in  Table  3,  and  PS~P/PS“1»  with  pg=1.88  Mg/m  . 

An  average  air  void  volume  of  7.0%  was  used  for  the  initial  air  void  content 
and  the  resulting  porous  crush  response  is  compared  in  Figure  31  to  that  of  the 
surrounding  tuff  and  the  crush  response  assumed  in  the  cylindrical  chimney  calculation. 
In  assuming  an  initial  void  content  of  7.0%,  the  bulking  effects  of  the  chimney 
collapse  process  have  been  neglected.  Core  data  and  reentry  observations  suggest 
that  this  approximation  may  be  quite  reasonable  in  the  lower  levels  of  the  chimney. 
Indeed,  more  recent  chimney  gas  flow  tests  suggest  that  the  average  air-filled 
porosity  of  the  chimney  as  a whole  lies  in  the  range  of  7.5  to  13  percent/20^  Thus 
any  tendency  for  increased  compaction  to  occur  at  lower  levels  in  the  chimney  must 
result  in  air-filled  void  volumes  at  or  near  the  assumed  7%  value.  These  average 
chimney  void  volume  measurements  are  also  consistent  with  the  drill  back  data^1^ 
which  suggest  the  conoidal  shape  assumed  in  the  calculations,  with  its  greatly 
increased  volume  in  comparison  to  the  cylindrical  chimney  geometry. 

As  indicated  in  Figure  30,  the  yield  surface  for  the  conoidal  chimney  material 
is  approximately  a factor  of  two  stronger  than  that  of  the  cylindrical  chimney  rubble 
at  the  same  confining  pressure.  Again  the  form  of  the  yie’d  surface  is  that  given 
by  Equation  (5)  with  constants  listed  in  Table  3.  Actually  the  difference  in  the 
effective  strengths  of  the  two  chimney  materials  is  much  greater  than  is  suggested 
by  Figure  30,  due  to  the  coupling  between  the  porous  crushup  and  yield  response 
models.  The  constitutive  model  used  to  represent  the  response  of  the  cylindrical 
chimney  material  is  extremely  dissipative  due  to  the  high  air  void  content  and  low 
pressure  required  for  void  collapse  (Figure  31).  As  a consequence,  in  this  material 
the  stress  induced  by  ground  shock  drops  very  rapidly  with  time.  Thus,  most  of  the 
motion  of  this  chimney  material  occurs  at  very  low  pressure  levels  where  the 
allowable  stress  difference  is  quite  small  (see  Figure  30).  In  the  much  more  competent 


Table  3 


Material  Properties  Used  to  Describe 
Tuff  and  Chimney  Rubble  Response  in  In-Chimney  Simulations 


Material 

Property 

Undisturbed 

Tuff 

Cylindrical 

Chimney 

Rubble 

Conoidal 

Chimney 

Rubble 

3 

Ambient  Density,  PQ  (Mg/m  ) 

1.920 

1.470 

1.750 

Poisson's  Ratio  , v 

0.33 

0.50 

0.33 

Bulk  Modulus,  Kq  (GPa) 

6.64 

0.214 

4.83 

Rigidity  Modulus,  Gq  (GPa) 

2.55 

0 to  1.04(a) 

1.85 

Compressional  Wave  Speed,  Cp  (m/s) 

2286 

382 

2043 

Shear  Wave  Speed,  Cg  (m/s) 

1152 

0 

1028 

Elastic  Pressure,  pg  (MPa) 

10 

10 

10 

Crush  Pressure,  pc  (MPa) 

400 

400 

370 

Air  Void  Content,  V (%) 

a 

2.0 

25.0 

7.0 

Unconfined  Yield,  Yq  (MPa) 

18.0 

0 

15.0 

Strength  Increase,  Ym  (MPa) 

52.0 

26.0 

39.0 

Maximum  Pressure,  p (MPa) 
m 

200 

120 

200 

r 

- 

- 

1.80 

A (GPa) 

- 

- 

5.0 

B (GPa) 

- 

- 

18.0 

(a)  Rigidity  modulus  varies  quadratically  with  distention  ratio,  a ; ranging 
from  0 to  1.04  GPa  as  pressure  increases  from  0 to  400  MPa,  respectively. 
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conoidal  chimney  material,  however,  stress  wave  attentuation  is  much  less  severe 
due  to  the  lower  void  content  of  this  material,  and  consequently,  the  effective 
strength  of  this  material  is  relatively  high  during  the  time  most  of  its  shock-induced 
motion  occurs. 

4.3  IN -CHIMNEY  STUDY  RESULTS 

4.3.1  Cylindrical  Chimney 

The  most  striking  feature  in  the  results  from  the  cylindrical  chimney  compu- 
tation involved  the  development  and  persistence  of  velocities  in  the  chimney  region. 
To  a lesser  extent,  the  focusing  which  is  produced  by  the  difference  in  material 
along  a "preferred  direction"  in  the  chimney  is  also  evident.  The  description  of 
some  of  the  field  variables  will  be  illustrated  through  the  use  of  isometric  plots. 
The  field  variables  are  plotted  up  vertically  from  the  spatial  plane  which  includes 
the  axis  of  symmetry  of  the  problem.  The  ranges  included  in  that  plane  will  be 
noted  on  the  figures,  as  will  the  scale  of  the  field  variable  axis.  The  reader  should 
note  that  since  these  plots  are  self-scaling,  the  amplitudes  on  the  field  variable  axis 
will  in  general  be  different  from  plot  to  plot  which  requires  some  "mental  scaling" 
to  be  performed.  In  considering  the  velocity  isometrics,  the  reader  should  note  that 
these  values  are  actually  the  speed  and  regions  of  negative  velocity  (towards  the 
detonation  point)  will  be  pointed  out  to  the  reader  where  they  are  significant. 

A comparison  of  the  velocity  fields  at  15  and  35  ms  (Figure  32)  shows  the 
shock  wave  developing  in  the  undisturbed  tuff  region,  as  well  as  in  the  chimney 
material.  The  difference  between  the  velocity  field  at  15  ms  and  that  at  the 
overlaid  time  of  almost  13  ms  is  some  additional  propagation  in  range  along  with  a 
drop  in  the  peak  velocity  from  455  meters  per  second  to  400  meters  per  second. 
By  35  ms,  it  is  clear  that  the  chimney  material  has  been  accelerated  to  significantly 
higher  velocity  than  the  undisturbed  material.  A cusp  in  the  wave  front  at  the 
cavity  boundary  Is  caused  by  the  relatively  slow  propagation  velocity  of  the  shock 
in  the  chimney  material.  By  50  ms  (Figure  33),  the  velocity  in  the  tuff  outside  the 
chimney  is  only  about  one-tenth  of  that  inside  the  chimney.  At  100  ms,  the  only 
significant  velocities  in  the  problem  are  those  in  the  chimney  material.  The  peak 
velocity  in  the  chimney  has  gone  down  by  more  than  a factor  of  2 from  that  at 


50  ms.  By  100  ms,  the  rebound  has  begun  and  the  region  next  to  the  cavity  in  the 
hemisphere  below  the  working  point  level  is  moving  back  towards  the  detonation 
point. 

Figure  34  shows  a profile  of  the  velocities  near  the  axis  of  symmetry,  also 
indicating  the  rebound  velocities.  A comparison  of  the  one-dimensional  spherical 
solution  conducted  entirely  in  chimney  material  with  the  two-dimensional  solution 
suggests  the  beginning  of  a focusing  effect.  In  the  one-dimensional  solution  since 
the  entire  sphere  is  expanding  at  the  same  rate,  the  cavity  pressure  drops  faster 
than  in  the  two-dimensional  problem  where  strong  tuff  lies  outside  the  chimney 
region  in  the  lower  hemisphere.  This  leads  to  a sustained  pressure  by  comparison 
which  continues  to  drive  the  chimney  material.  It  is  already  clear  that  the  residual 
stress  fields  formed  at  later  times  will  be  grossly  affected  in  this  region  since 
chimney  material  is  not  rebounding  with  the  surrounding  tuff. 

For  an  explosion  of  this  yield  in  the  undisturbed  medium,  the  motions  would 
ordinarily  be  minimal  at  about  190  ms.  After  that  time,  relatively  small  oscillations 
around  zero  endure  in  velocity.  The  half  period  for  such  oscillations  in  this  calculation 
is  about  60  ms.  However,  the  chimney  region  material  continues  to  move  outward 
from  the  detonation  point  as  shown  in  Figure  35.  At  200  ms,  the  motion  in  the 
free-field  is  predominantly  back  towards  the  detonation  point,  but  the  motion  in  the 
chimney  region  is  outward.  Owing  to  the  faster  propagation  velocity  in  the  surrounding 
tuff,  the  wave,  by  coupling  through  to  the  region  at  the  top  of  the  chimney,  has 
actually  produced  motions  of  chimney  material  back  towards  the  cavity  (Figure  36). 
It  is  also  seen  in  the  figure  that  a large  region  of  the  chimney  continues  to  move 
outward  at  250  ms  with  a maximum  velocity  of  about  20  meters  per  second. 

By  300  ms,  the  material  at  the  top  of  the  chimney  moving  inwards  has  stacked 
up  against  the  enormous  mass  of  material  still  moving  away  from  the  detonation 
point.  The  peak  velocity,  which  is  relatively  uniform  in  the  chimney  region,  is  now 
close  to  9 meters  per  second.  The  chimney  material  continues  to  slow  down  with 
the  majority  of  it  moving  at  about  3 meters  per  second,  but  with  a peak  velocity 
of  6 meters  per  second  near  its  leading  edge  (Figure  37).  By  350  ms,  the  material 
is  starting  to  meet  resistance  from  the  surrounding  tuff  at  the  top  of  the  chimney 
region.  The  computation  was  ended  at  this  time.  A more  quantitative  illustration 
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Isometric  speed  plot  at  250  ms  for  cylindrical  chimney. 


of  the  outward  velocity  in  the  chimney  at  350  ms  is  given  in  Figure  38.  Notice 
that  the  steep  part  of  the  wave  has  just  reached  the  top  of  the  chimney  at  a range 
of  150  meters.  It  is  clear  that  the  stagnation  of  the  chimney  material  at  the  top 
of  the  chimney  will  produce  high  stresses  there  which  will  propagate  back  toward 
the  cavity.  The  extent  of  these  "water-hammer"  like  oscillations  will  depend  on  the 
ability  of  the  undisturbed  material  adjacent  to  the  chimney  to  constrain  the  lateral 
motion,  as  well  as  the  boundary  condition  provided  by  the  cavity  pressure. 

The  residual  stress  field  in  the  vicinity  of  the  chimney  is  still  changing  at 
350  ms  when  the  computation  was  terminated.  However,  in  the  region  adjacent  to 
the  detonation-produced  cavity,  the  residual  stresses  have  eventually  established 
themselves  by  200  ms.  The  evolution  of  the  hoop  stress  component  from  200  ms 
to  350  ms  is  illustrated  in  the  contour  plots  shown  in  Figure  39.  A region  of  stress 
below  that  of  the  original  overburden  of  70  bars  is  seen  to  grow  in  extent  even 
producing  tensile  regions  within  it.  This  region,  adjacent  to  the  chimney,  would 
ordinarily  be  of  low  residual  stress  and  it  appears  that  the  relief  into  the  chimney 
region  has  further  reduced  that  stress  producing  a valley  through  the  ordinarily  high 
residual  stress  field  around  the  cavity.  This  valley  seems  to  be  directed  back  toward 
the  detonation  point  level.  This  would  clearly  unacceptable  for  the  safety  for  HLOS 
tests  in  which  a tunnel  would  ordinarily  be  placed  at  working  point  level.  The  valley 
in  stress  would  be  a possible  path  of  communication  around  the  region  which  is 
ordinarily  stemmed.  An  isometric  plot  of  the  hoop  stress  at  250  ms  is  shown  in 
Figure  40.  This  shows  that  the  magnitude  of  the  average  hoop  stress  in  the  chimney 
region  is  below  that  of  the  cavity  pressure.  The  valley  in  the  hoop  stress  referred 
to  above  is  seen  to  be  well  below  the  chimney  hoop  stress  and  directed  back  towards 
the  detonation  point  elevation. 

The  reduction  of  stress  due  to  the  presence  of  the  chimney  is  shown  by  the 
profiles  in  Figure  41.  These  profiles  are  taken  along  lines  perpendicular  to  the  axis 
of  symmetry  at  different  distances  above  and  below  the  detonation  point  level.  The 
hoop  stress  profiles  olotted  below  the  working  point  level  are  very  similar  to  those 
computed  for  the  pancake  problems  discussed  in  Section  3 (see  Figure  24)  when  the 
factor  of  two  in  yield  is  considered.  At  only  12  meters  above  the  detonation  point 
level,  a 30  percent  reduction  is  noted  in  the  peak  stress. 
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Figure  38.  Axial  velocity  near  the  axis  of  symmetry  at  350  ms  for  cylindrical 
chimney. 
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4.3.2  Conoidal  Chimney  and  Comparison 

The  major  differences  between  the  two  in-chimney  computations  concern  the 
chimney  geometries  and  the  properties  of  the  chimney  material  (Section  4.2).  In 
terms  of  conducting  parameter  studies,  this  is  somewhat  unfortunate  since  in  general, 
it  is  difficult  to  identify  differences  caused  by  the  various  changes.  However,  after 
reviewing  the  results  it  became  clear  that  the  substantive  differences  in  the  compu- 
tations were  caused  by  the  differences  in  material  properties,  as  opposed  to  those 
in  geometry.  These  differences  are  discussed  in  what  follows. 

The  conoidal  chimney  problem  was  initialized  from  the  one-dimensional  compu- 
tation at  7.69  ms.  It  will  be  remembered  that  the  cylindrical  chimney  problem  was 
initialized  at  a time  of  12.84  ms,  the  difference  being  due  to  the  propagation  velocity 
in  the  two  different  chimney  materials.  By  50  ms,  the  ground  shock  has  established 
itself  well  out  into  the  undisturbed  tuff,  as  well  as  up  into  the  chimney  material 
(Figure  42).  By  comparison  with  Figure  33,  it  is  seen  that  the  velocity  surface  for 
the  conoidal  chimney  is  significantly  more  symmetric  about  the  detonation  point 
than  was  the  case  for  the  cylindrical  chimney.  The  cylindrical  chimney  showed  peak 
velocities  above  the  cavity  which  were  more  than  a factor  of  2 greater  than  those 
below  the  cavity,  and  also  on  the  axis  of  symmetry. 

At  a time  of  100  ms  when  the  rebound  in  the  region  adjacent  to  the  lower 
hemisphere  of  the  cavity  is  establishing  itself,  the  chimney  material  has  been 
accelerated  to  a peak  velocity  of  about  17  meters  per  second  (Figure  43).  The  form 
of  the  velocity  profile  on  axis  is  similar  to  the  cylindrical  chimney  case;  however, 
the  peak  velocity  for  the  cylindrical  chimney  was  more  than  four  times  greater  at 
80  meters  per  second.  By  200  ms,  the  highest  velocity  in  the  conoidal  chimney 
case  is  almost  an  order  of  magnitude  down  from  that  of  the  cylindrical  case  (compare 
Figure  44  with  Figure  35).  The  shape  of  the  profile  in  the  conoidal  chimney  also 
exhibits  a dip  at  the  center  at  that  region  which  is  significantly  less  well  developed 
in  the  cylindrical  chimney  case.  Some  additional  wave  structure  is  evident  for  the 
conidal  chimney  (Figure  44).  Some  of  this  appears  to  be  real  and  has  been  caused 
by  the  greater  coupling  between  the  chimney  material  and  the  undisturbed  material 
for  this  case.  The  structure  is  more  difficult  tc  see  on  Figure  35  since  the  scaling 
for  the  isometric  plot  has  been  chosen  to  accommodate  the  very  high  peak  velocity 
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A 


A. 


Isometric  speed  plot  at  50  ms  for  conoidal  chimney  problem. 


Isometric  speed  plot  at  100  ms  for  conoidal  chimney  problem.  Material  below  cavity  is 
rebounding,  while  chimney  material  continues  to  move  outward. 


rebounding  material  exist. 


on  axis.  Some  elements  of  chimney  material  are  already  moving  back  towards  the 
cavity  in  the  conoidal  chimney  case  at  this  time.  This  effect  was  also  noted  for 
the  cylindrical  chimney. 

By  250  ms,  the  magnitude  of  the  peak  velocities  is  almost  exactly  an  order 
of  magnitude  lower  for  the  conoidal  chimney;  however,  this  velocity  is  in  a direction 
back  towards  the  cavity.  The  peak  velocities  from  the  cylindrical  case  were  still 
moving  outward  at  about  20  meters  per  second  over  a substantial  length  of  the 
chimney  with  only  a small  region  at  about  2 meters  per  second  moving  back  towards 
the  cavity.  For  the  conoidal  chimney  case,  these  low  velocities  continue  to  rattle 
around,  but  are  no  longer  playing  a role  in  establishing  the  residual  stress  field  which 
is  discussed  below. 

The  most  striking  difference  between  the  two  chimney  computations  is  that 
for  the  conoidal  chimney,  the  residual  stress  fields  which  were  established  at  250 
ms  are  not  significantly  different  from  what  one  would  expect  for  a computation 
without  a special  chimney  region.  The  residual  hoop  stress  is  shown  at  a time  of 
250  ms  in  Figure  45  for  the  conoidal  chimney  case.  Note  the  near-symmetric 
magnitude  of  the  hoop  stress  peaks  around  the  cavity.  It  is  seen  that  the  hoop 
stress  in  the  chimney  region  above  the  cavity  is  greater  than  that  in  surroundings 
and  that  the  gradient  moving  away  from  the  detonation  point  is  somewhat  lower 
than  the  region  beneath  the  cavity.  At  250  ms  for  the  cylindrical  chimney  case, 
the  residual  stresses  in  the  direction  above  the  cavity  were  essentially  completely 
eroded  contrary  to  what  is  seen  here  for  the  conoidal  case.  This  is  what  made  the 
tensile  "valley"  in  the  cylindrical  case  appear  to  be  ominous.  Although  such  a valley 
also  exists  for  the  conoidal  case,  the  fact  that  the  residual  stresses  around  the 
cavity  are  so  high  casts  the  evaluation  of  the  low  stress  region  in  a substantially 
different  light. 

A more  quantitative  illustration  of  the  hoop  stress  in  different  regions  at  250 
ms  for  the  conoidal  problem  is  shown  in  the  contour  plot  of  Figure  46.  Notice  the 
particular  asymmetry  of  the  region  bounded  by  the  40  MPa  contour.  This  region, 
which  is  of  finite  extent  beneath  the  cavity,  tapers  away  to  a cusp  before  reaching 
the  axis  of  symmetry  above  the  cavity.  Another  asymmetry  is  that  of  the  10  MPa 
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Isometric  plot  of  residual  hoop  stress  (at  250  ms)  for  conoidal  chimney  problem. 


contour  which  shows  a kink  at  the  chimney  boundary.  This  is  the  same  type  of 
effect  which  was  seen  on  the  cylindrical  chimney  problem,  but  is  of  significantly 
less  importance  here.  The  small  tensile  region  adjacent  to  the  chimney  also  appears 
more  clearly  in  the  contour  plot.  In  the  cylindrical  calculation  this  region  was  the 
beginning  of  the  low  stress  valley  to  which  reference  was  made.  In  the  conoidal 
case,  the  cavity  is  clearly  isolated  from  the  tensile  region  by  the  development  of 
the  residual  stress  field  which  has  formed. 

4.3.3  Conclusions 

The  outstanding  conclusion  from  the  in-chimney  study  is  that  the  properties 
of  the  material  which  fill  the  chimney  region  are  of  paramount  importance  regarding 
the  formation  of  a competent  residual  stress  field.  The  materials  chosen  here  in 
one  sense  represent  two  bounds  of  what  is  likely  to  exist  in  an  actual  chimney. 
The  distribution  of  cavity  volume  throughout  the  chimney  in  the  cylindrical  case 
provides  an  upper  bound  on  what  the  air  void  content  of  the  recompacted  material 
at  a chimney  bottom  is  likely  to  be.  While  this  may  underestimate  the  amount  of 
void  in  the  upper  levels  of  the  chimney,  it  clearly  overestimates  the  amount  of  void 
in  the  lower  levels.  The  use  of  properties  from  overlying  strata  is  probably  reasonable 
in  the  cavity  region,  but  underestimates  the  air  voids  in  the  rubble  at  higher  levels 
of  the  chimney. 

Because  of  the  clear  dependence  of  the  residual  stress  field  on  the  material 
properties  for  in-chimney  configurations,  one  would  have  to  carefully  document  the 
chimney  properties  during  the  exploratory  phase  of  site  selection.  Extensive  work 
to  document  these  properties  would  have  to  be  conducted  on  a case  by  case  basis 
in  order  to  qualify  any  particular  site.  Before  the  actual  exploration,  it  would  be 
prudent  to  perform  computations  which  included  the  effects  of  layered  media  both 
in  the  chimney  and  the  surrounding  tuff. 
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